M = 512; % 原图像长度
N = 64; % 水印图像长度
K = 32; % 子块大小

alpha=0.05;% 嵌入强度系数

% 打开原图、水印图
I = imread('Cameraman.bmp');
G = imread('lsc_sign.jpg');
%W = zeros(M);
 
% 缩放、灰度化原图、改变精度
I = imresize(I,[M M]);
%I = im2double(I); % double精度转换
% I = rgb2gray(I); % 灰度化处理
 
G = imresize(G,[N N]);
%G = im2double(G); % double精度转换
G = rgb2gray(G); % 灰度化处理
 

subplot(5,4,1);
imshow(I);
title('原始载体图片');
subplot(5,4,2);
imshow(G);
title('原始水印图像');
 

%Step 1
[LL,LH,HL,HH] = dwt2(G,'haar'); % 进行2维haar离散小波变换
[U,S,V] = svd(HH);% 对HH进行SVD分解，得到U、S、V矩阵
 
%Step 2
%进行2级离散小波变换
[LL1, LH1, HL1, HH1] = dwt2(I, 'haar');
[LL2, LH2, HL2, HH2] = dwt2(LL1, 'haar');%128*128
H0 = entropy(HH2);% 计算HH3系数的信息熵
 
%Step 3
%选出最优嵌入块 默认为4*4:(1,1)
optimal_block_index = 0;

%Step 4
%对最优嵌入块进行 DCT 变换，得到DCT系数矩阵 B
m = floor(optimal_block_index/4)+1;
n = mod(optimal_block_index, 4)+1;
x = (m - 1) * K + 1;
y = (n - 1) * K + 1;
H_I = HH2(x:x+K-1, y:y+K-1);
B = dct2(H_I);








% %Step 5
% %对B进行奇异值分解,嵌入水印
% [U1,S1,V1] = svd(B);
% S2 = S1 + alpha * S;
% B1 = U1 * S2 * V1;
% H_I = idct2(B1);
% HH2(x:x+K-1, y:y+K-1) = H_I;
% LL1 = idwt2(LL2,LH2,HL2,HH2,'haar');
% W = idwt2(LL1,LH1,HL1,HH1,'haar');
% W = uint8(W);

%Step 5 带arnold变换
% 对B进行奇异值分解,嵌入水印
[U1, S1, V1] = svd(B);
S2 = S1 + alpha * S;
B1 = U1 * S2 * V1;
H_I = idct2(B1);
HH2(x:x+K-1, y:y+K-1) = arnold(H_I, 1, 1, 5); % 使用Arnold变换更新HH2
LL1 = idwt2(LL2, LH2, HL2, HH2, 'haar');
W = idwt2(LL1, LH1, HL1, HH1, 'haar');
W = uint8(W);


% 计算载水印图像和宿主图像之间的PSNR，评估水印算法不可感知性
psnr_img = psnr(I, W);
fprintf('PSNR about img: %f.\n', psnr_img);


%攻击
%高斯低通滤波攻击 
W_GAOS = gaos(W, 0.5);
%JEPG 压缩攻击
W_JEPG = jepg(W);
% 剪切攻击
W_CROP = crop(W, 0.8);
% 旋转攻击
W_ROTA = rotate(W, 13);
% 高斯噪声
W_GSN = gaus_noise(W, 0.5);
% 椒盐噪声
W_SP = sp(W, 0.5);
% 泊松噪声
W_POSI = posi(W);
% 乘性噪声
W_X = xxing(W, 0.5);

% 提取水印 原图，高斯低通，压缩，裁剪，旋转，高斯噪声，椒盐噪声，泊松噪声，乘性噪声
A = get_watermark(W, x, y, K, S1, U, V, LL, LH, HL, alpha);
A_GAOS = get_watermark(W_GAOS, x, y, K, S1, U, V, LL, LH, HL, alpha);
A_JEPG = get_watermark(W_JEPG, x, y, K, S1, U, V, LL, LH, HL, alpha);
A_CROP = get_watermark(W_CROP, x, y, K, S1, U, V, LL, LH, HL, alpha);
A_ROTA = get_watermark(W_ROTA, x, y, K, S1, U, V, LL, LH, HL, alpha);
A_GSN = get_watermark(W_GSN, x, y, K, S1, U, V, LL, LH, HL, alpha);
A_SP = get_watermark(W_SP, x, y, K, S1, U, V, LL, LH, HL, alpha);
A_POSI = get_watermark(W_POSI, x, y, K, S1, U, V, LL, LH, HL, alpha);
A_X = get_watermark(W_X, x, y, K, S1, U, V, LL, LH, HL, alpha);

subplot(5,4,3);
imshow(W);
title('嵌入水印后的载体图像');
subplot(5,4,4);
imshow(A);
title('提取出来的水印图像');

subplot(5,4,5);
imshow(W_GAOS);
title('高斯低通滤波攻击后的载体图像');
subplot(5,4,6);
imshow(A_GAOS);
title('高斯低通滤波攻击后提取的水印');

subplot(5,4,7);
imshow(W_JEPG);
title('JEPG压缩攻击后的载体图像');
subplot(5,4,8);
imshow(A_JEPG);
title('JEPG压缩攻击后提取的水印');

subplot(5,4,9);
imshow(W_CROP);
title('随机剪切攻击后的载体图像');
subplot(5,4,10);
imshow(A_CROP);
title('随机剪切攻击后提取的水印');

subplot(5,4,11);
imshow(W_ROTA);
title('旋转攻击后的载体图像');
subplot(5,4,12);
imshow(A_ROTA);
title('旋转攻击后提取的水印');

subplot(5,4,13);
imshow(W_GSN);
title('高斯噪声攻击后的载体图像');
subplot(5,4,14);
imshow(A_GSN);
title('高斯噪声攻击后提取的水印');

subplot(5,4,15);
imshow(W_SP);
title('椒盐噪声攻击后的载体图像');
subplot(5,4,16);
imshow(A_SP);
title('椒盐噪声攻击后提取的水印');

subplot(5,4,17);
imshow(W_POSI);
title('泊松噪声攻击后的载体图像');
subplot(5,4,18);
imshow(A_POSI);
title('泊松噪声攻击后提取的水印');

subplot(5,4,19);
imshow(W_X);
title('乘性噪声攻击后的载体图像');
subplot(5,4,20);
imshow(A_X);
title('乘性噪声攻击后提取的水印');


% 计算PSNR值 所有攻击提取出的水印的PSNR
% fprintf('PSNR about watermark: %f.\n', psnr(G, A));
% fprintf('PSNR about GAOS: %f.\n', psnr(G, A_GAOS));
% fprintf('PSNR about JEPG: %f.\n', psnr(G, A_JEPG));
% fprintf('PSNR about CROP: %f.\n', psnr(G, A_CROP));
% fprintf('PSNR about ROTATION: %f.\n', psnr(G, A_ROTA));
% fprintf('PSNR about GAOS_NOISE: %f.\n', psnr(G, A_GSN));
% fprintf('PSNR about SP: %f.\n', psnr(G, A_SP));
% fprintf('PSNR about POSITION: %f.\n', psnr(G, A_POSI));
% fprintf('PSNR about Xxing: %f.\n', psnr(G, A_X));




% 计算NC值：体现提取受攻击后的水印质量
nc = caculate_nc(G, A);
fprintf('NC about watermark: %f.\n', nc);
nc_gaos = caculate_nc(G, A_GAOS);
fprintf('NC about GAOS: %f.\n', nc_gaos);
nc_jepg = caculate_nc(G, A_JEPG);
fprintf('NC about JEPG: %f.\n', nc_jepg);
nc_crop = caculate_nc(G, A_CROP);
fprintf('NC about CROP: %f.\n', nc_crop);
nc_rota = caculate_nc(G, A_ROTA);
fprintf('NC about ROTATION: %f.\n', nc_rota);
nc_gsn = caculate_nc(G, A_GSN);
fprintf('NC about GAOS_NOISE: %f.\n', nc_gsn);
nc_sp = caculate_nc(G, A_SP);
fprintf('NC about SP: %f.\n', nc_sp);
nc_posi = caculate_nc(G, A_POSI);
fprintf('NC about POSITION: %f.\n', nc_posi);
nc_x = caculate_nc(G, A_X);
fprintf('NC about Xxing: %f.\n', nc_x);


% 下面是函数————————————————————————————————————————————
% NC计算函数
function nc_val = caculate_nc(G, A)
% 计算直方图
h1 = imhist(G);
h2 = imhist(A);
% 根据直方图计算 NC 值
nc_val = sum(sqrt(h1 .* h2)) / sqrt(sum(h1) * sum(h2));
end



% 提取水印 带arnold变换
function A = get_watermark(W, x, y, K, S1, U, V, LL, LH, HL, alpha)
[LL3, ~, ~, ~] = dwt2(W, 'haar');
[~, ~, ~, HH4] = dwt2(LL3, 'haar'); % 128*128
H_I2 = HH4(x:x+K-1, y:y+K-1);
B2 = dct2(H_I2);
[~, Sw, ~] = svd(B2);
Sx = (Sw - S1) / alpha;
B2 = U * Sx * V;
H_I2 = idct2(B2);
A = idwt2(LL, LH, HL, rearnold(H_I2, 1, 1, 5), 'haar'); % 使用Arnold变换还原
A = uint8(A);
end

%高斯低通滤波攻击
function W_GAOS = gaos(W, sd)
H = fspecial('gaussian',3, sd);
W_GAOS = imfilter(W,H);
end
%JEPG 压缩攻击
function W_JEPG = jepg(W)
quality = 50;
W_JEPG = imresize(W, 0.5); % 缩小图像
imwrite(W_JEPG, 'temp.jpg', 'Quality', quality); % 保存为JPEG格式
W_JEPG = imread('temp.jpg'); % 重新读取JPEG图像
W_JEPG = imresize(W_JEPG, 2); % 放大图像
end
% 剪切攻击
function W_CROP = crop(W, r)
% r = 0.3; % 剪切比例为30%
sz = size(W);
W_CROP = W;
h1 = round(sz(1)*r); % 剪切高度
w1 = round(sz(2)*r); % 剪切宽度
x1 = round(rand(sz(1)-h1)); % 随机选择一行
y1 = round(rand(sz(2)-w1)); % 随机选择一列
W_CROP(x1+1:x1+h1, y1+1:y1+w1) = 0; % 将指定区域置为0
end
% 旋转攻击
function W_ROTATE = rotate(W, angle)
% angle = 20; % 旋转角度为20度
W_ROTATE = W;
W_ROTATE = imrotate(W_ROTATE, angle, 'bilinear', 'crop');
end
% 高斯噪声
function W_GSNOISE = gaus_noise(W, r)
W_GSNOISE = imnoise(W, 'gaussian', 0, r);
end
% 椒盐噪声
function W_SP = sp(W, r)
W_SP = imnoise(W, 'salt & pepper', r);
end
% 泊松噪声
function W_POSI = posi(W)
W_POSI = imnoise(W, 'poisson');
end
% 乘性噪声
function W_X = xxing(W, r)
W_X = imnoise(W, 'speckle', r);
end
% arnold变换代码
%img 灰度图像 a,b为参数 n为变换次数
function arnoldImg = arnold(img,a,b,n)
[h,w] = size(img);
N=h;
arnoldImg = zeros(h,w);
for i=1:n
    for y=1:h
        for x=1:w
            %防止取余过程中出现错误，先把坐标系变换成从0 到 N-1
            xx=mod((x-1)+b*(y-1),N)+1;
            yy=mod(a*(x-1)+(a*b+1)*(y-1),N)+1;  
            arnoldImg(yy,xx)=img(y,x);              
        end
    end
    img=arnoldImg;
end
arnoldImg = uint8(arnoldImg);
end

% 逆arnold变换代码
function img = rearnold(arnoldImg,a,b,n)
[h,w] = size(arnoldImg);
img = zeros(h,w);
N = h;
for i=1:n
    for y=1:h
        for x=1:w           
            xx=mod((a*b+1)*(x-1)-b*(y-1),N)+1;
            yy=mod(-a*(x-1)+(y-1),N)+1  ;      
            img(yy,xx)=arnoldImg(y,x);              
        end
    end
    arnoldImg=img;
end
img = uint8(img);
end